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ABSTRACT 

# A multi-stage Runge— Kutta method is analyzed for solving the two- 

dimensional Euler equations for external and internal flow problems. 

Subsonic, supersonic and, highly supersonic flows are studied. Various 

techniques for accelerating the convergence to a steady state are described 
and analyzed. Effects of the grid aspect ratio on the rate of convergence are 
evaluated. An enthalpy damping technique applicable to supersonic flows is 
described in detail. Numerical results for supersonic flows containing both 
oblique and normal shocks are presented confirming the efficiency of the 

method. 
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I. INTRODUCTION 


Over the past few years, a number of numerical schemes have been devised 
for solving the compressible Euler equations. An explicit finite-volume scheme 
with multi-stage Runge-Kutta integration in time, as proposed by Jameson, 
Schmidt, and Turkel [1], has shown promise of being an efficient and accurate 
Euler solver. This method has been found to be of great versatility in the 
simulation of transonic flow, particularly when the method has been augmented 
by such convergence accelerating techniques as local time-stepping, enthalpy 
damping, and implicit residual averaging. The incorporation of these 
techniques as described by Turkel [2] has produced remarkable increases in the 
rate of convergence to a steady-state. However, the extent to which these 
enhancements are necessary vary significantly from one problem to another. 
Although the preceding method has been used successfully for transonic flow, 
its use in high-speed flows involving shocks of substantial strength has not 
been explored fully. Moitra [3] has presented some results for supersonic 
flow demonstrating the possibility of its use in capturing complex shock 
structures. The aim of this study is to investigate the applicability, 
accuracy, and efficacy of the multi-stage Runge-Kutta method in solving the 
compressible Euler equations for two-dimensional supersonic internal flow 
problems involving complex shock discontinuities of substantial strength. The 
qualitative and quantitative nature of the effects of various techniques for 
accelerating the convergence to a steady-state are analyzed in detail. The 
results from the present calculations are compared with experimental data and 
results obtained by other numerical methods [4,5]. 



II. SOLUTION ALGORITHM 


The Runge-Kutta time-stepping scheme is more efficient than many other 
explicit schemes such as MacCormack's original scheme. In addition, the 
steady-state solutions are independent of the time-step, and as mentioned 
earlier, the scheme is amenable to a variety of techniques for accelerating 
steady-state convergence. Before the results are presented, we present a 
description of these accelerating techniques. 

Local Time-Stepping 

A variable time-step determined by the bound on the local Courant number 
is used In each zone. The use of the largest time-step at any location allowed 
by the local stability bound accelerates convergence to the steady-state at 
the expense of time-accuracy of the solution. 


Implicit Residual Averaging 

This strategy for accelerating convergence consists of taking an average 
of residuals in neighboring cells by the following expression 


eRj_j + (l-2e)R i + eR i+i ~ R i 


( 1 ) 


where denotes the residual in the i-th cell and is the residual 

before the implicit averaging. This is done after every even stage of the 
scheme. This results in a scheme which is stable for a CFL number X if 


e 




1 


( 2 ) 
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A 

where X is the CFL limit of the original scheme. Since X can be much 

higher than X , the rate of convergence is accelerated. However, choosing 

too large a CFL limit, X, slows down the convergence. An optimal choice is 
& 

X/X ~ 2. See reference [2] for further details. 

Artificial Dissipation 

The method requires artificial dissipative terras in order to alleviate 
two difficulties associated with centrally differenced schemes. The first one 
is that of odd-even point decoupling that adversely affects convergence. 
Addition of a fourth-difference artificial viscosity reduces this difficulty. 
The second difficulty manifests itself in oscillations in the vicinity of 
shocks. An additional viscous term is introduced, which is an approximation to 
a second derivative in order to overcome this problem. To maintain second- 
order accuracy this terra contains a coefficient that itself depends on the 
second difference of the pressure. Since the fourth-order difference can 
create oscillations near the shock, it is turned off in the vicinity of 
shocks • 

Enthalpy Damping 

The enthalpy damping procedure has been modified from the original 
version [1] of FL052 in order to enable it to be used in supersonic flow. 
Hence, we shall describe the new enthalpy damping in more detail. 

In solving the full potential equation it is well known that adding 
artificial terms such as <)> t can accelerate the convergence to a steady 
state. For a first-order system of equations such as the Euler equations, 
this can not be done. Instead, it is suggested [1] that one can add forcing 
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functions that depend on h-tig, h = (E+P)/p; where h and hp denote the 

local and the free stream values of the total specific enthalpy. In the 

steady state, the total specific enthalpy h is constant along each 
streamline in an inviscid flow. Assuming that the flow comes from a common 
reservoir, in the steady state h = hQ everywhere, where Viq is determined 
by the inflow boundary conditions. Based on the analogy with the full 
potential equations and the fact that h — hp « <f> , it was suggested that the 
forcing function be added to the density equation. Using density p and 
velocity components (u,v) the modified equations are 


P t + L i = -ap(h - h Q ) 


(2a) 


u t + l 2 = 0 


(2b) 


v t + l 3 = °f 


(2c) 


where Lj are the usual space derivatives of the Euler equations. When one 
solves for the momentum variables pu and pv, then the forcing function 
also occurs in these equations so that 


(pu) t + L 2m = " apu ( h ' V 


(3a) 


(pv) t + L 3m = _apv ( h " V 


(3b) 


where the subscript m denotes the space derivative terms in the momentum 
equations. The finite difference approximation to the density equation is 



p n = -At - aAt p n+1 (h n - h Q ) 


( 4 ) 


n+1 

P 


or 


n+1 

P 


P n - At Lj 
1 + aAt(h n - h Q ) 


(5) 


Freezing coefficients and changing variables [2], one can reduce (2) to a 
second-order equation for p, 


P TT - (c 2 - u 2 )p nn + 2Kp x = 0 (6) 

where n and t are new independent variables and the positive constant K 

n 

depends on 1 - M* and linearly on a. Thus, we obtain the potential-like 
equation with the artificial term Kp . To see the effect of K, we assume 
that p has the simple form 


P = e® T -7r < n < 7r. (7) 

Assuming subsonic flow, letting q 2 = c 2 - u 2 and substituting (7) into (6), 
we find 

2 2 2 
6 + 2KB + C » 0 

or 

6 = -K + / K 2 - q 2 £ 2 . (8) 

In order to reach a steady state as fast as possible we wish the real part of 
3, i.e, Re 3 to be as negative as possible. Note for K = 0, i.e., no 



enthalpy damping, Re (3 " 0 and we have no term driving the system to a steady 
state. We see from (8) that Re 3 is most negative when 0 < K < q£ , and 
K is as large as possible. Thus, we wish to choose K slightly less than 

q?. 

We first consider the case in which we want to damp the low frequencies, 
since these are the slowest to decay. In this case, we have 5 = 0(1) and 
so K = 0(1). We stress that the constant HM in the FL052 code is equal to 
aAt. Hence, if a = 0(1), this implies that HM should vary linearly with At 
and therefore, is very small. We next consider the case where we wish to damp 
the high frequencies, e.g., multigrid applications. Now £ = tt/Ati and so 
K varies as 1/An or equivalently 1/At. In this case, HM = aAt should be 
a constant. 

This analysis implies that the enthalpy damping suggested in [1] is 
useful if: 

(1) the flow is locally subsonic; 

(2a) HM is proportional to At to damp low frequencies; 

(2b) HM is constant to damp high frequencies. 

In order to explain the effectiveness of the enthalpy damping for non- 
raultigrid applications and to extend its usefulness to supersonic applications 
we must discuss the energy equation which was ignored in (2). Since the 
enthalpy damping was originally based on an analogy with the full potential 
equation there is no obvious way to treat the energy equation. The original 
idea was to append to (2) the equation 


ds 

at 


+ 


J 4s 


- 0 , 


(9) 
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i.e., nothing special is done to the entropy equation. Rewriting (9) for the 
energy equation we obtain 


If + L 4E " h o) h ‘ 

However, now the forcing function contains a quadratic term in h and was 
found to be destabilizing rather than stabilizing. To avoid this difficulty 
an ad hoc correction was used and the energy equation became 

If + L 4E = " a ( h " h o)‘ (2d) 

Instead of the ad hoc procedure, we shall replace (9) by 

H +L 4h = - 6 ( h - h o)' <9 ' ) 

Thus, instead of ignoring the entropy we try and drive the total enthalpy to 
its steady state value. We note that if 3 = 00 we recover the isoenergetic 
approximation in which the energy equation is replaced by Bernoulli's law. 
Thus, our new system is a compromise between the original Euler equations and 
the isoenergetic equations. Of course, in the steady state all of these 
variations are equivalent. Eliminating h we can rewrite (9') as 

||+ L 4E = - (ctE +|£)(h - h Q ). (2d') 

We thus, conclude that enthalpy damping helps in the non-multigrid version of 
FL052 because of the energy equation and not because of the continuity 
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equation, i.e., we can set a = 0 in (2) and the coefficient 3 assumes the 
important role. This improved enthalpy damping can be used equally well in 
supersonic regions. Furthermore, one can choose 3 as large as one desires 

without causing instabilities, while choosing too large an aAt causes 

instabilities. However, choosing 3 too large may slow down the convergence 
rate. In all the cases presented in the results section we chose a = 0 and 
3At to be constant. 

Combining all the previous discussions we prescribe for the system (2a), 
(3a), (3b), and (2d") that: 

(1) For locally supersonic flow we must have a * 0. For locally subsonic 

flow we advise choosing a to depend on 1 - M^. The coefficient 3 in 

(2d") can be used for both subsonic and supersonic flows and does not depend 
on the Mach number. 

(2) To dampen the low frequencies, HM = aAt should depend linearly on 

At while to damp the high frequencies, HM = aAt should be chosen to be 
constant. The choice of 3 is frequency independent, but BAt should be 
constant. Hence, even for a multigrid version one should choose aAt and 
3At independently, especially since a should depend on 1 - M . 


III. RESULTS 

A wide variety of test cases involving different Mach numbers and shock 
strengths were investigated in order to assess the ability of the numerical 
scheme to simulate flows ranging from subsonic to high supersonic Mach 
numbers. Convergence histories for the different test cases are presented 
with a view to establishing the nature and extent of the effects of the 
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different convergence accelerating techniques as well as different grid 
sizes. The CPU-seconds execution times presented are those for scalar 
processing on a CDC Cyber-175 processor. Computed results have been compared 
with experimental data and results from other schemes whenever available. 

Test Case 1 

The first case investigated is the inviscid shock reflection problem in 
which an incident shock wave of prescribed strength reflects off a flat 
plate. A schematic of the configuration is presented in Figure 1. Uniform 
grids with various numbers of points in the x-direction (0 < x < 4.0) and 
the y-direction (0 < y < 1.5) were considered. At the inflow boundary for 
y > 1 and along the upper boundary all dependent variables are specified 
based on a prescribed incident shock-angle a. All variables are extrapolated 
from the interior at the supersonic outflow boundary. 

The first set of computations for this configuration was performed with a 
grid of 51 points in the x-direction and 31 points in the y-direction. The 
shock-incidence angle a was 29° • A contour plot of the computed pressures 
appears in Figure 2. Both the incident and the reflected shocks are seen to 
be reasonably sharp. In Figure 3, a comparison of the pressure profile at the 
y " 0.25 line computed by the present method with results obtained by the 
original MacCormack scheme [4] and the Osher [5] upwind scheme is presented 
and a good agreement of the results is noticed. Next, a set of convergece 
histories corresponding to the use of the different convergence accelerating 
techniques, is presented in Figure 4. The use of a locally maximum time step 
is seen to substantially increase the rate of convergence over that obtained 
by using a global minimum time step. The use of enthalpy damping does not 
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create any advantage for this particular configuration; however, the 
beneficial aspects of enthalpy damping will be apparent in the results 
presented for the next test case. Finally, the use of the implicit residual 
averaging technique results in a slight further increase in the convergence 
rate. In all these plots the logarithm of the residual, i.e., the average 
value of the time derivative of density, is plotted against actual CPU seconds 
of computation. 

The next series of computations for this configuration were done in order 
to ascertain the effect of the grid aspect ratios on the rate of convergence. 
Seven uniform grids of different sizes were used to compute the flow with 

= 2.9 and a = 29° and the correspondinormack scheme [4] and the Osher 
[5] upwind scheme is presented and a good agreement of the results is 
noticed. Next, a set of convergece histories corresponding to the use of the 
different convergence accelerating techniques, is presented in Figure 4. The 
use of a locally maximum time step is seen to substantially increase the rate 
of convergence over that obtained by using a global minimum time step. The 
use of enthalpy damping does not create any advantage for this particular 
configuration; however, the beneficial aspects of enthalpy damping will be 
apparent in the results presented for the next test case. Finally, the use of 
the implicit residual averaging technique results in a slight further increase 
in the convergence rate. In all these plots the logarithm of the residual, 
i.e., the average value of the time derivative of density, is plotted against 
actual CPU seconds of computation# 

The next series of computations for this configuration were done in order 
to ascertain the effect of the grid aspect ratios on the rate of convergence. 
Seven uniform grids of different sizes were used to compute the flow with 
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= 2.9 and a = 29° and the corresponding convergence histories are 
presented in Figure 5. The different aspect ratios of the cells are seen to 
have pronounced effects on the convergence rates. The fastest rates of 

convergence were generated by three particular grids, namely those with 
41x25, 51x31, and 61x37 points. These three grids have a cell aspect ratio 
of approximately 1.62. A departure from this value of the cell aspect ratio 
is seen to substantially slow the rate of convergence. This indicates that 
the convergence of the finite-volume Euler scheme is strongly influenced by 
the aspect ratio of the computational cells. We next consider the possibility 
that the angle of the shock may have some influence on the value of the cell 
aspect ratio for which the fastest rate of convergence is attained. To test 
this hypothesis an additional series of computations were performed for the 
same configuration with the shock angle a changed to 23°, all other 
parameters remaining the same. The resulting convergence histories are 
plotted in Figure 6. As is immediately obvious, the grid sizes with cell 
aspect ratio of 1.62 again give rise to the fastest convergence rates. A 
departure from this value, as in the case of the 51x43 grid, is again seen 
to result in substantial worsening of the rate of convergence. This leads to 
the conclusion that moderate changes in the shock angle do not substantially 
influence the value of the optimal cell aspect ratio for which the fastest 
convergence rate is achieved. 

Test Case 2 

The next test cases involved internal flow in a single throat nozzle [6]. 
Flow regimes ranging from subsonic to highly supersonic in nature have been 
investigated. A variety of complex shock-structures have been successfully 
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captured. At supersonic inflow boundaries all variables are specified and at 
supersonic outflow all the variables are extrapolated from the interior. For 
subsonic inflow boundaries characteristic-type boundary conditions are 
applied. At outflow boundaries characteristic extrapolation is supplemented 
by a radiation boundary condition. A uniform gird of 61x31 points is used 
in all of these computations. 

A fully supersonic internal flow is considered as the first test problem 
with this configuration. In this problem Mach 3.5 flow enters the two- 
dimensional nozzle which is symmetric about the centerline. A contour plot of 
the computed pressures is presented in Figure 7. A complex shock structure 
consisting of multiple— ref lected shocks is noticed. Interaction of the 
expansion fan emanating from the nozzle wall with the shock reflected off the 
centerline is clearly visible in the throat region. This demonstrates the 
ability of the present Euler method to capture repeatedly reflected shock 
structures at moderately high supersonic Mach numbers. The effects of the 
various convergence accelerating strategies are apparent in the convergence 
histories plotted in Figure 8. Enthalpy damping, in this case substantially 
increases the rate of convergence. The use of implicit residual averaging is 
seen to remarkably accelerate convergence and causes an approximately 200 
seconds reduction in the computational time required to reach the same level 
of the residual as that in the case of enthalpy damping alone. 

The next problem provides a more severe test of the solution procedure's 
a bility to capture shocks and to handle subsonic and supersonic flow regions 
in the configuration. The incoming flow is subsonic at a Mach number of 
0.46. The flow accelerates as it travels down the nozzle and becomes super- 
sonic at the throat creating a supersonic region in the middle section of the 
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nozzle. This supersonic region ultimately terminates at a normal shock beyond 
which the flow is again subsonic. The terminal shock is clearly visible in 
the pressure— contour plot presented in Figure 9. The presence of these 
different flow regimes do not necessitate any special treatments as far as the 
algorithm is concerned. A comparative plot of surface pressure computed by 
the present algorithm and available experimental data [6] is presented in 
Figure 10. The agreement of the results is noticeably good and confirms the 
ability of the method to predict complex internal flows accurately. Some 
experiments at hypersonic Mach numbers have been performed by the authors and 
the results obtained generated encouragement as to the applicability of the 
method to computation of hypersonic flow. 


IV. CONCLUSIONS 

Enthalpy damping, heretofore restricted in its applicability to subsonic 
flows, has been demonstrated to be successfully applicable to supersonic 
flows. Significant savings in computational time have been achieved by the 
use of the present enthalpy damping technique. Problems involving supersonic 
and subsonic flow regions have been solved successfully demonstrating the 
versatility of the Runge-Kutta scheme. A particular value of the cell aspect 
ratio is seen to cause the fastest rate of convergence; however, moderate 
changes of the shock angle are not seen to have significant influence on this 
value of the aspect ratio. 
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Figure 3. Plot of computed pressure at y = 0.25 for 2-D shock 
reflection. M = 2.9, a = 29°. 
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Figure 4. Convergence histories for different accelerating techniques for 2-D 
shock-reflection. = 2.9, a = 29°, Grid: 51 x 31. 
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Figure 5. Convergence histories for different grid sizes for 2-D shock 
reflection. M = 2.9, a = 29°. 

* denotes cell aspect ratio - 1.62. 
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Figure 7. Contour plot of computed pressure for 2-D nozzle. M^ n = 3.5. 
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Figure 8. Convergence histories for different accelerating techniques for 
2-D nozzle* M. = 3*5. 
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Figure 10. Comparison of computed and experimental surface pressure for 2-D 

nozzle. M. - 0.46. 
in 





1. Report No. NASA CR 178053 ^ overnment Accession No. 

ICASE Report No. 86-5 

3. Recipient’s Catalog No. 

4. Title and Subtitle 

APPLICATION OF A RUNGE-KUTTA SCHEME 
FOR HIGH-SPEED INVISCID INTERNAL FLOWS 

5. Report Date 

.T^nn^ry 1986 

6. Performing Organization Code 

7. Author(s) 

Anutosh Moitra, Eli Turkel, and Ajay Kumar 

8. Performing Organization Report No. 

86--5 

9. Performing Organization Name and Address 

Institute for Computer Applications in Science 
and Engineering 

Mail Stop 132C, NASA Langley Research Center 
4Jamoton VA ^3665‘“522 t ^ 

10. Work Unit No. 

11. Contract or Grant No. 

NAS1-17070; NAS1-18107 

13. Type of Report and Period Covered 
Contractor Report 

12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washi npTnn . P . C . 20546 

14. Sponsoring Agency Cone 

505-31-8^-01 

U ' — ^ — — ' 

15. Supplementary Notes 

Langley Technical Monitor: AIAA Paper No. 86-0104 

J. C. South AIAA 24th Aerospace Sciences 

Meeting, 1/6-9/86, Reno, NV 

JFirml Paprrrt-_ 

16. Abstract 


A multi-stage Runge-Kutta method is analyzed for solving the two- 
dimensional Euler equations for external and internal flow problems* 
Subsonic, supersonic and, highly supersonic flows are studied* Various 
techniques for accelerating the convergence to a steady state are described 
and analyzed. Effects of the grid aspect ratio on the rate of convergence are 
evaluated. An enthalpy damping technique applicable to supersonic flows is 
described in detail. Numerical results for supersonic flows containing both 
oblique and normal shocks are presented confirming the efficiency of the 
method. 


17. Key Words (Suggested by Authors(s)) 


18. Distribution Statement 


Euler, 


34 - Fluid Mechanics & Heat Transfer 

Runge-Kutta, 
high-speed, 
internal flow 


Unclassified 

- unlimited 


19. Security Classif.(of this report) 

20. Security Classif.(of this page) 

21. No. of Pages 

22. Price 

Unclassified 

Unclassified 

26 

A0 3 


For sale by the National Technical Information Service, Springfield, Virginia 22161 


NASA-Langley, 1986 




















